A many-body approach to crystal field theory 
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A self-consistent many-body approach is proposed to build a first-principles crystal field theory, 
where crystal field parameters are calculated ab initio. Many-body theory is used to write the 
energy of the interacting system as a function of the density matrix of the noninteracting system. 
A variation of the energy with respect to the density matrix gives an effective Hamiltonian matrix 
that is diagonalized to determine the density matrix providing the lowest energy. The equations are 
written in terms of the quantum group of functional derivatives with respect to external fermionic 
sources. This approach contains the many-body theory of Green functions as a special case and the 
usual crystal field theory as a first approximation. Therefore, it is expected to provide good results 
for strongly-interacting electron systems. 



The many-body theory of realistic quantum systems 
uses two main methods: (i) the Green function method 
that gives good results for extended states such as metals 
and semiconductors and (ii) the diagonahzation method, 
typically used in quantum chemistry, which describes ef- 
ficiently localized systems and the splitting of degenerate 
states. The latter model is efficient to calculate the mag- 
netism of rare earths and transition metal impurities in 
insulators, as well as of their optical and spectroscopic 
properties. However, the second method leads to the 
diagonahzation of very large matrices for complex sys- 
tems. Thus, an intermediate approach is desirable, where 
a small-size Hamiltonian is set up using effective param- 
eters. The best-known example of this approach is the 
crystal-field or ligand-field method, which was extended 
by Kotani and coll. to describe the interaction be- 
tween a localized and an extended system. It gives good 
results for the calculation of x-ray absorption and pho- 
toemission spectroscopies. 

The main drawback of this method is its use of parame- 
ters, which can be numerous in the case of low symmetry. 
The calculation of crystal field and Slater parameters in 
molecules and solids is a subject of continued interest 
H H H H 0, 0, II ES lll^, but they are not part of 
a general formalism and double counting of interactions 
is difficult to control. Moreover, the effect of the crystal 
field is not taken into account self-consistently. 

In this paper, we propose a first-principles self- 
consistent approach of the effective Hamiltonian method, 
based on the quantum field theory of correlated systems 

[Ulli. 

It was noticed a long time ago by Esterling and Lange 
that the many-body solution of a degenerate or 
quasidegenerate system contains a loophole: among the 
various (quasi) degenerate states of the noninteracting 
system, only one will evolve towards the ground state of 
the interacting system, and it must be known in advance 
to set up the many-body calculation. Here we solve this 
problem by showing that the relevant state of the non- 
interacting system can be obtained by diagonalizing an 



effective Hamiltonian. Moreover, as shown in \\'^ - the re- 
quirement of self-consistency imposes the use of density 
matrices instead of pure states. So, for a density matrix 
p, we calculate the total energy of the interacting system 
E{p). To determine the density matrix p that minimizes 
the energy, we have to diagonalize an effective matrix 
which is determined by the quantum group approach to 
many-body theory. Therefore, this method unifies the 
diagonahzation method and the Green function method. 

Total energy of the system - Standard many-body the- 
ory 01 enables us to write the total energy of an interact- 
ing system which evolves from a noninteracting system 
described by the density matrix p as 

Eip) = tr(p5-iT(if(0)e-*''-^==^"(*)^*)), (1) 

where S is the S-matrix, T is the time-ordering operator 
and H{t) = Ho{t) + ff™'(i) is the total Hamiffonian in 
the interaction representation, with 



Hoit) 



2 Vni,r)^t(i,i.')^(i,r')V(t,r) 



drdr', 



i/-*(t) = - 

where /io(r) = —A /2m + {/^{r) is the one-body Hamil- 
tonian operator for an electron in a nuclear potential 
C/jv(r). 

Within the functional derivative approach to quantum 
field theory 0, this expression can be rewritten 
E{p) = PiZp), where 



P(Zp) = flimhoir) 



(5?7+(x)(5r/+(i/) 



dr 



-drdr'. 



2|r — r'l Sr]+ {y)6ri+ {x)Sr]+ {x)6r]+ (y) 

where x = (i,r), y = {t,r') and the right-hand side is 
taken for zero external sources. Zp is the generating func- 
tion of the Green functions. It is given by the expression 
Zp — J2kl PlkZkl, with plk the density matrix and 



= {K\S{r,rr'S{fi+,^+)\L). 
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Here, S{r]^,r]^) is the S-matrix in the presence of 
fermionic external sources fj^ and 77^ and Zkl — 
e-'-°Z^^, where D = D+ - D_ and D± is 
the interaction term ^^^H™^{t)dt where the fields 
are replaced by functional derivatives with respect 
to the external sources fj±,ri±. Finally, Z^j^ — 
exp[— i J f]{x)GQ{x,y)r]{y)dxdy]N^j^. In this equation, 
Tj is a two-dimensional vector with components 77+ and 
r]-, Go(a:,y) is the 2x2 matrix 



z(0|V(x)V^(y)|0) 



-i{0\^l;Hy)^j{x)\0) 
0\T*{^{x)^l;Hy))\0) 



T* is the anti-time-ordering operator and iV]^^ — 
{K\:exp[i J f]d{x)il){x) + {x)rid{x)Ax]:\L) , with rjd = 
77+ - 77_. 

To calculate N^j^, we must define the states \K) and 
\L) . The solution of the noninteracting Schrodinger equa- 
tion provides one-electron orbitals Un(r) with energy e„. 
An orbital is called a core orbital if it is filled in all states 
\K) and |L), otherwise, it is called a valence orbital. The 
core orbitals are numbered from 1 to C, the valence or- 
bitals from 1 to M. There are C electrons in the core 
orbitals and N < M electrons in the valence orbitals. 
For example, in the ion Cr'^"'", there are C = 18 core or- 
bitals and = 3 electrons in the M — 10 orbitals of the 
degenerate 3d shell. The states are generated from the 
vacuum |0) by the action of creation operators cjj for core 
electrons and vl for valence electrons as 



\K) 







= • 





■cIlO), 
.cl|0), 



where ik, jk are valence orbitals (i.e. integers taken in 
the set {1, . . . , M}) ordered so that ii < ■ ■ ■ < ipf and 
ji < • • • < jAT. A lengthy calculation 13] yields 



= ]^(l + afeQ;fe) 



k=l 



M 



exp(^ 



92 



dundar, 



with an = J rid{x)un(x)dx and q;„ — J u'^ nix)rid{x)dx, 



(r) and 



where u„(x) 
X = {t,r). 

Notice that, for vanishing sources (i.e. when 7^^ — 
ff^ — 0), Zkl — Sk.l- To calculate E(p), it is also 
possible to use the Galitskii-Migdal formula [i2| , but our 
expression for E(p) allows for a more direct comparison 
with the standard equations of the crystal field theory. 

Energy minimization - The minimization of E{p) with 
respect to p is constrained by the fact that p must be 
Hermitian, nonnegative and with unit trace. To remove 
the trace condition, we use the fact that, for a gen- 
eral matrix p, Zp = tr(p) for vanishing sources. So 



the trace condition is relaxed by writing the energy as 
E{p) = P{Zp)/Zp. The conditions of Hermiticity and 
nonnegativeness are relaxed by writing p — b^b, where 
b is now unconstrained. If we take the derivative of 
E{p) — P{Zp)/Zp with respect to 6^ and look for an 
extremum where dE{p) / dh^ = wc obtain the equa- 
tion E{p)dZp/db^ = P{dZp/db^) or, more exphcitly 
E{p)Y.m^kmZml = Y.m^kmP{Zml)- For vanishing 
sources we have Zml = 5l,m and this equation becomes 
E{p)b = bP{Z), where Z is the matrix Zml- 

To solve this, we diagonalize P{Z). This gives us a 
unitary matrix U and a diagonal matrix z such that 
E{p)b = bUzU\ or E(p)bU = bUz. 

To determine the general solution of this equation, we 
order the eigenvalues and eigenvectors of P{Z) and we 
call eo the smallest eigenvalue. We put R — bU, and we 
assume that the lowest eigenvalue is d-fold degenerate. 
The general solution of E{p)bU = bUz (i.e. E{p)R = Rz) 
giving the minimum energy is E{p) — eo and Rij = for 
j > d. Then b = RW and p ^ UR^RW, where R is an 
arbitrary Hermitian nonnegative d-dimensional matrix. 
At this stage, we can determine R by considering the 
symmetry of the system. The solution which gives the 
same weight to each eigenstate (and then the maximum 
symmetry) is R = 1^. This determines the density 
matrix to be 



PKL 



d 

E 



(2) 



where the factor 1/d was added to have tr(p) = 1. The 
maximum symmetry can be broken by choosing another 
d-dimensional matrix R^ R. For example, in a spheri- 
cally symmetric system, the matrix R^ R can be chosen 
to select a specific spectroscopic term "^^^^L. In sys- 
tems where the spin-orbit interaction is weak, we can 
select a given spin state. This can be useful to make self- 
consistent calculations of low-spin and high-spin transi- 
tion metal ions. 

Tie effective matrix - It remains now to determine the 
effective matrix P{Z). We can obtain P{Z) by noticing 
that Zp = J2klPlkZkl implies Zkl = dZp/dpLK- 
Therefore, the equation for Zkl can be obtained from 
an equation for Zp by a derivation with respect to plk- 

The equation for Zp is [l3| 



dZ, 



(3) 



where 9 is a functional derivative with respect to an ex- 
ternal source 77^ (x) or 77^(0;). In fact, we need an infinite 
system of equations (the Green function hierarchy) ob- 
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tained by taking the functional derivative of equation 

ddZp = ^^^^(-l)-^(dili?il'"aWpi)(di2Di2'"Zp), 

(4) 

where / = |Dil'"| + |di2||Dil'"| + |di2|, d is a polynomial 
in the functional derivatives with respect to an exter- 
nal source f]^{x) or ?]^{x) and all quantities are taken 
for vanishing sources (i.e. fj — rj ^ 0). In the rest of 
the discussion, we consider only the equation for dZp, 
the other equations being obtained through further func- 
tional derivatives, as for equation 

Taking the partial derivative of equation ||2J) with re- 
spect to plk, we obtain an equation for OZkl- 

^^^^ - 5]^$](-l)l^'i'"l((I?il'"5W^^J(i5i2'" 

n— 

+ (i^il'"aTypi)pi2'"ZKL), 

where W]^i^ — dWp / dpLK- For the calculation of the 
partial derivative of Wp with respect to plk, it is neces- 
sary that a general value of p be used in Wp . In particu- 
lar, p must not be assumed Hermitian or with unit trace. 
This impHes the presence of terms such as l/iv{p) in Wp. 
Once the partial derivative is calculated, the condition 
tr(p) = 1 can be restored to simpHfy W}^]^. Notice that 
Wp is not linear in p, so M^^^ depends on p. This depen- 
dence implies that the solution of the system of equations 
for p, Zp and Zkl must be solved self-consistently. 

We are interested in the energy splitting due to the 
interaction. So we can remove the terms that are pro- 
portional to 5kl because they just shift all eigenvalues 
by a constant term and do not modify the eigenvectors. 
To do that we write Zkl = Ykl + {Zp + Y)5kl^ with the 
boundary conditions Ykl = and Y — for vanishing 
external sources. Introducing this into equation 10), we 
obtain the following equation for Ykl'- 

^^^^ - f;^5](-l)l^'^'"l(pil'"9yKL)(i?i2'"Zp) 

+ {Dil"'dWl){Di2"'YKL)), (6) 

where 

VkL = WkL-P'T.^kK, 

dim 

with dim the dimension of the matrix W]^j^. 

Finally, the algorithm for the calculation of the density 
matrix p and the Green functions dZp is the following: (i) 
choose an initial density matrix p, (ii) calculate Wj, and 
Vk l , (iii) use equation to calculate the dZp you need 
for the calculation of the energy, (iv) calculate d Yr- l for 
a given value of dZp and Vkl, (v) diagonalize P{Ykl), 



(vi) deduce a new density matrix Pnew, (vii) compare 
Pnow and p, if convergence is not reached, go back to (ii) 
with a new p. Remark that steps (iii) and (iv) require 
the solution of an infinite system of equations. This sys- 
tem must be broken by some closure approximation of 
the Green function hierarchy and the resulting dZp and 
dYKL are thus approximate. Notice also that another 
self-consistency loop can be introduced in the algorithm 
to calculate a self-consistent external potential, as was 
described in This will change the orbitals Un{x). 

The crystal field equations - We must show that 
the present approach contains the standard crystal field 
equations as a first-order approximation. The func- 
tions we need to calculate P{Ykl) are OO'Ykl with 
d' — S/6ri^{y) and d = S/Si]^{x) and did2dzdiYKL 
with di = 5/5r]+{y), 83 = S/Sr]+{x), 82 = S/6fj+{x) 
and di — S/6fj~^{y). If we solve iteratively equations Q 
^p) and lO for 71 = we obtain 

(5) dd'YKL - dd'VKL, (7) 

did2d3d4YKL did2d3d4VKL - {d2d4VKL){dld3Zp) 
+ {dld4VKL){d2d3Zp) - {d2d4W^){did3YKL) 
+ {dld4W^)(d2d3YKL)- (8) 

A lengthy calculation shows that P{Ykl) calculated from 
(O and ^ gives the same matrix as the standard crystal 
field theory. 

Tize simplest example - As an example, we consider the 
case of one electron (N=l) that can occupy two orbitals 
(M=2). The density matrix is then of dimension 2. The 
matrix Vkl takes the form 

_ f (aiai - 612012)12 0201 \ 

\ 0x02 [a.2012 — aiQ!i)/2 J 

( {Pii - 912)1^ P12 \ 

V P2I (P22-Pll)/2y 

In particular, if the system is degenerate (pn = P22 = 
1/2, P12 — P21 — 0), the second term is zero and Vkl has 
only a term linear in a and o. 

We presented a first-principles self-consistent many- 
body approach to crystal field theory. Our approach 
contains by construction the many-body Green function 
theory as the case N = M and it contains the crystal 
field theory as a first order approximation. Thus, it is 
expected to provide good results in situations that are 
intermediate between the range of validity of these two 
limits, such as in strongly-interacting electron systems. 

This method can be used more generally to calcu- 
late the splitting, due to interaction, of levels that are 
degenerate in the noninteracting system, for any quan- 
tum field theory. A interesting example is the applica- 
tion of this approach to the quantum electrodynamics of 
atoms, where its self-consistency should make it compet- 
itive with respect to the effective Hamiltonian methods 
proposed by Shabaev and Lindgren [l9| . 
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I am very grateful to S. Di Matteo for mentioning 
the important reference Q|. This is IPGP contribution 
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